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APPLICATION OF A FINITE DIFFERENCE TECHNIQUE 
TO THERMAL WAVE PROPAGATION 

by Kenneth J. Baumeister 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

ABSTRACT 

A finite difference formulation is presented for thermal wave propa- 
gation resulting from periodic heat sources. The numerical technique 
can handle complex problems that might result from variable thermal 
diffusivity, such as heat flow in the earth with ice and snow layers. 

In the numerical analysis » the continuous temperature field is repre- 
sented by a series of grid points at which the temperature is separated 
into real and imaginary terms. Next, computer routines previously de- 
veloped for acoustic wave propagation are utilized in the solution for 
the temperatures. The calculation procedure is illustrated for the 
case of thermal wave propagation in a uniform property semi- infinite 
medium. 

INTRODUCTION 

Because of increasing concern over environmental heat losses, en- 
vironmental heat balances are now an integral part of many design pro- 
cedures. In certain applications, environmental heat transfer calcu- 
lations can be concerned with the effect of periodic solar temperature 
variations. An example of such a problem might be the impact of oil 
pipe lines on the periodic temperature distribution beneath the arctic 
tundra. In order to analyze problems of this type as well as more 



2 


general thermal wave propagation problems, a finite difference formu- 
lation is developed in this paper. 

At present, analytical solutions of thermal wave propagation 
(refs. 1 and 2) are limited to relatively simple geometries with con- 
stant thermal properties or lumped constant temperature regions. The 
present finite difference technique, however, can handle the heat flow 
field complications which might arise from (1) variations in thermal 
properties, (2) complex geometries, or (3) complex boundary conditions. 
For example, the flow of heat into the earth with covering ice and 
snow layers could be readily calculated. 

The finite difference technique presented herein represents a 
direct extension of a difference theory developed for noise propagation 
in ducts (refs. 3 and 4). First, the temperature field is represented 
by a series of grid points in which the temperature at each grid point 
is separated into real and imaginary terms. Next, the governing energy 
equation and the appropriate boundary conditions are present and 
written in difference form. Finally, computer routines previously 
developed for acoustic wave propagation are utilized for the solution 
of the resulting thermal wave difference equations. The calculational 
procedure is illustrated for the case of thermal wave propagation in 
a uniform property semi-infinite medium. 
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SYMBOLS 

A coefficient matrix, eq. (34) 

a thermal diffusivity 

C submatrix, eq. (35) 

c specific heat 

F matrix eq^ (34) 

i V”1 

k thermal conductivity 

m total number of grid rows 

n total number of grid columns 

T temperature 

T q amplitude of temperature boundary condition 

temperature in far field 
t time 

X axial position of exit plane 

x dimensionless axial coordinate, eq. (9) 

Ax axial grid spacing 

Y transverse width of temperature field 

y dimensionless transverse coordinate, eq. (10) 

Ay transverse grid spacing 

Z exit impedance, eq. (23) 

0 dimensionless temperature 

0^ boundary condition 

p density 
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spatial temperature field 
<t> o boundary condition 

to circular frequency 

Superscripts : 

1 dimensional quantity 

J 1 or 2 

K 2 or 1, K=J- (-1) J 

(1) real part 

(2) imaginary part 
Subscripts : 

i, j i, axial index, j, transverse index, see Figure 2 

r reference quantity 

1 material //l 


2 


material #2 
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GOVERNING EQUATION AND BOUNDARY CONDITIONS 
Conduction Equations 

The governing differential equation describing the conduction of 
heat energy in its two dimensional form is given by 


9 9T __8_ , 8T _ 8T 

8x T 8x T 8y ’ 8y T PC 9t 


( 1 ) 


where the prime, 1 , is used to denote a dimensional quantity. Define 
a dimensionless temperature 

T - T 

e = T --r - f • " . 

o 00 

where T^ is the initial temperature of the medium assumed to be uniform 
and T q is the amplitude of the temperature variation at the surface 
or other appropriate convenient temperature. Therefore, the governing 
equation (1) becomes 


8 , 99 89 __ _86 

8x T *8x f 3y ’ k 8y f pC 3t 


(3) 


Equation (3) can be rewritten as 


1_ __8_ _80_ 1_ _J_ 89 = 1__ 89 

pca^ 3x* k 8x' pca^ 8y' 8y T 8t 


where the thermal diffusivity is defined as 


a pc 


( 4 ) 


( 5 ) 
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and the reference diffusivity 


a = 


r pc 
r r 


( 6 ) 


For a periodic heat source (neglecting the initial transient) , 
the solution for the temperature can he assumed to be of the form 


(x 1 ,y 1 ,t) = <J>(x,y)e 


iwt 


(7) 


Substituting equations (5) and (7) into equation (4) gives 


a 3 


- k _ii + i_ k J*£ _ ± m 


ka^ 3x 9x' ka^ 9y f 3y' 


<j> = 0 


( 8 ) 


Finally, the spatial coordinates x 1 and y' are non-^dimensionalized 
by assuming 


x' = 


2tt 


r \l 2 

— 1 
2a J 

v r / 


(9) 



( 10 ) 
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Thus , 


+ i 87T 2 (f, = 0 

ka^ dx 3x ka^ 3y 3y r • 


( 11 ) 


For uniform properties (a equals a^) , aquation (II) reduces to 


+ Btt 2 * = 0 (12) 

3x 3y 

This equation is independent of the frequency parameter co. The fre*- 
quency parameter enters the problem by coordinate stretching the 
dimensionless variables x and y. 

Using the exponential notation displayed in equation (7), the di^ 
mensionless temperature $ has , in general, both real and imaginary 
parts . 

Thus , 


(J»(x,y) = 4> (1) + i^ (2) 


(13) 


Consequently, equation (12) can be broken into Its real and imaginary 
parts by substituting equation (13) into equation (12) , 


2,(1) 2 .(1) m 
■5-A— - + -Ll_ + 87I 2 ^ 2 ) = o 

3x 3y 


, 2 . { 2 ) „2 ( 2 ) 
±±-+*-± T -- Sir 2 / 1 * = 0 
3x Z 3y 


(14) 


( 15 ) 



The above two equations can be expressed as a single equation 


aV J) , a¥ J) 

3x 2 3y 2 


(-l) J 8ir 2 <j> (K) 


= 0 


( 16 ) 


where J equals 1 or 2 and K equals 2 or 1 respectively. 

The solution for temperatures can be found by either solving the 
complex equation (12) or the set of real equations represented by 
equation (16). From the standpoint of computer efficiency, the com- 
plex equations are the easiest to handle. Equation (16) will only 
be used to indicate what type of solution schemes are most advantageous. 
Therefore, any remaining equations or boundary conditions will be 
left in complex form. Of course, the true physical temperature is 
represented by only the component. 

Boundary Conditions 

Surface temperature . - The surface temperature (see Figure 1) 
treated in this paper assumes 


0(O,y,t) = 0 Q (y) cos ait 
which in complex form is written as 

0(0, y,t) = 0 o (y)e 1Wt 

Therefore, the surface condition on $ becomes 


(17) 


( 17 ) 
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where $ is a real quantity. 

Adiabatic . - The temperature distribution at y = 0 and Y is 
chosen such that the heat flux is zero. Thus, 


3y 


o 



(19) 


Exit temperature . - The exit condition chosen sets the mean level 
about which the temperature oscillation can occur. For example, the 
exit temperature could be assumed to be of the form 

T « at x = X (20) 


0 = 0 at x = X 


( 21 ) 


Exit impedance . - Another possible exit condition is very useful 
in treating a problem in which the medium is assumed to be semi- infinite 
in the x direction. Consider the one dimensional problem in which 
the exit plane is at infinity. The analytical solution is derived in 
reference 1 (pg. 329) and shown to be (for uniform properties) 


-2tt(1+1)x 

e 


( 22 ) 


Borrowing an analogy from acoustics, the ratio between the temper- 
ature (J) and the temperature gradient is defined as the impedance 

Z exit ; 
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Z exit 


3x 


For a one- dimensional semi-infinite medium, 


Z 


exit 



(23) 


(24) 


The ratio defined by equation (24) is not a function of x. Thus, the 
temperature field of a semi-infinite medium can be closely approximated 
by a finite medium, if equation (24) is used as the exit boundary 
condition. Of course, for a one-dimensional semi-infinite problem, the 
application of equation (24) at any position will give an exact agree- 
ment (within numerical accuracy) with analysis, as will be shown later 
in an example. 

SOLUTION TECHNIQUE 
Finite Difference Equations 

Instead of a continuous solution for temperature, the temperature 
will be determined at isolated grid points by means of the finite 
difference approximations, as shown in Figure 2. This enables changing 
the differential equations to a system of algebraic equations for the 
temperature at each grid point. The governing equations can be approxi- 
mated in difference form (ref. 5) by using either a Taylor series 
expansion, a variational, or an integral formulation. In this problem, 
where a gradient is specified along the exit boundary, the integral 
method for generating the finite difference approximation is most convenient. 
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The energy equation in finite difference form is developed by 
applying the integration method (ref. 5, pg. 168) to the cells in 
Figure 3. The cells are enclosed by the dashed lines which are 
spaced midway between the grid lines (not shown) . The grid lines would 
go directly through the grid points. Thus, the integration of the 
energy equation, equation (12), over the cell becomes 



r* —i 

2 „2 

rr i 

+ ^4 - i 

‘-"cell j 

9x 9y Z 


dxdy 


(25) 


By applying Green’s theorem for the plane region, equation (25) becomes 



dxdy = 0 


(26) 


where dS is a length element along the unit cell boundary, and N is 
the outward normal to the cell boundary. 

As usual, the difference formulation approximates the first 
derivative 


90 

9x 





90 



*1*3+1 ^i»3 

Ay 


(27) 


( 28 ) 


and so forth. 
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Cell 1 - central . - The integration equation, (26), is now 
applied to the central region of the temperature field, labeled cell 1 
in Figure 3. In evaluating the surface integrals in equation (26), 
the gradients are assumed to be constant on each side of the cell. 
Thus, equation (26) becomes 


2 


/ \2 

/ Ax\ 6. 1 

— l, j-1 

v A y / ■ 

+ <!>.-, . - 2 

/Ax\ 2 

1 + + i (2 ttAx) 

[AyJ 



+ .4> 


i+l,j 


+ 



i,j+l 


= 0 


(29) 


Some intermediate algebraic steps involved in the derivation of equation 
(29) and the following equations are presented in Appendixes C and D 
of reference 3, for a difference equation very similar to equation (29), 
Cell 2 - adiabatic boundary . - Now consider the difference -equation 
which applies in Cell 2, which is adjacent to the upper boundary in 
figure 3. For this unit cell, equation (26) can be expressed as 



,m-l 



* • -i 

i-l ,m 


1 + 



+ i(2?rAx) 2 


i,m 


+ I *1+1, m " 0 < 30 > 

where equation (19) was used to evaluate the term along the adiabatic 
boundary. A similar equation applies at the lower boundary. 

Cell 3 - exit plane . - In a similar manner, the difference equation 
which applies to cell 3 is found to be 
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1 

2 



n,j 


•i + Vi.j 


1 + 



+ i (2-rrAx) 2 + (l+i)2TrAx 





n,j+l 


0 


(31) 


where equation (24) written in the form 


3 4> 

9x 


n,j 


4tt 

i-1 9 n,j 


(32) 


was used to evaluate at the exit. 

3N 

Cell b ~ corner . - Finally, the difference equation for cell 4 is 


A a , A 4 A 

2 \Ayi ™n,m-l 2 ui-l,m 2 


1 + I™) + 1 (2ttAx) 2 


+ (1+i) 2 it Ax 


<f> - 0 

n,m 


(33) 


MATRIX SOLUTION 

The collection of the various difference equations at each grid 
point forms a set of simultaneous equations which can be expressed in 
matrix notation as 


{A} - [$] = [F] (34) 

where A is the known coefficient matrix, lj> is the unknown column vector 
containing the unknown complex temperatures, and [F] is the known column 
vector containing the various initial conditions. 



14 


Equation (34) can also be expressed in terms of all real 
quantities . In order to accomplish this , the column vector is 
expressed in terms of and and the A matrix is subdivided as 

follows 


r 


p - 


- - 

i 

i i 

n 

i 

\ 

i ( «_ 


jd) 

]+C j A 

l 1 

v J 

! 

* (2) 


~(2)~ 


The has a form typical of those matrices found in two-dimensional 
steady state heat-conduction problems, while the C matrix represents 
the coupling that occurs between the $ ^ and temperatures in the 

(-1)^ term of equation (16). 

Standard computer texts such as reference 5 dealing with the steady 
state heat conduction problem point out that the matrix is positive 
definite in at least one row; consequently, iteration schemes such as 
Gauss-Seidel can be used to obtain a solution. However, in the matrix 
equation (35), the C matrix adds to the off-diagonal elements ai*d as a 
result the matrix will no longer be positive definite. As a result, 
conventional iteration techniques cannot be used; However, matrices 
of the form of equation (34) or (35) can be solved by elimination 
techniques. In particular, a solution of the block tridiagonal form of 
the complex matrix equation (34) appears to be the most efficient from 
both operational time and computer storage consideration. Dr. D. W, Qqinn 
of Wright-Patterson Air Force Base, Dayton, Ohio has developed a 
computer package which will efficiently handle equation (35). This 
code had to be modified slightly to give results In terms of temperature 
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rather than acoustic pressure for which it was developed. 

APPLICATIONS 
Semi-Infinite Body 

As an example of the use of the numerical technique, the case of 
one-dimensional heat wave propagation into a semi-infinite medium will 
now be considered. This case allows a comparison of the numerical and 
analytical temperature distributions in the medium. The analytical 
solution for this problAm is given by equation (22) in this report. 

Because there is no variation of temperature in the y direction, 
the two-dimensional grid lattice shown in Figures 2 and 3 can be 
reduced to a one-dimensional lattice as shown in the upper sketch of 
Figure 4. The calculation was made for a. depth of 1 with the exit 
impedance given by equation (24) to simulate a semi-infinite medium. 

As seen in Figure 4, the agreement between the numerical and 
analytical results is good. By a series of numerical calculations, 
the number of grid points necessary to get accurate . temperature 
variations for this example was found to be 

n~12X (36) 

where X represents the position of the exit plane as shown in Figure 1, 

Composite Structure 

Finally, in order to solve problems with varying composition, the 
material can be assumed to be broken up into many bands each with a 
uniform composition for which equation (29) applies. However, to 
complete the problem, it is necessary to develop the difference equation 
at the junction between two materials as shown in Figure 5. 
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integration of equation (11) over the cell shown in Figure 5 gives 

+ +, 



— + r 5 - I- - i Sit 2 


ka^ 9x 3x ka^ 3y 3y 


dxdy - 0 


(37) 


The integration limits - and + stands for the left and right sides (also 
bottom and top) of the cell respectively. Integration of equation (37) 
gives 



(38) 


It is necessary, as shown above, to break the x integration into two 
domains (- to 0 and 0 + to +) , since the integrand is discontinuous 
at the center of the cell, labeled 0 in Figure 5. Now, equation (38) 
can be directly integrated to give 



2 

- i 87 t tfiy AxAy = 0 


(39) 



17 


which yields 


ifeiVMf * + !i* 

2 V a r / l fly / a r 


£* * 6=4 fe 


+ i 2(2 ttAx) 


5 i»j + \ K j ^ 1+1,3 


/ a. +a«\ a \ 2 
. 1 1 2 /Ax\ _ 

+ 2 l a i ( Ay ) ° 


(40) 


where it has been assumed that 


3tj) 
l X 3x 


- a 


8 (ft 
2 9x 


0 + 


(41) 


The approximation sign is used because the true interface boundary 
condition is that the heat flux are equal at the interface f Equation (39) 
can now be used a$ the interface between two different materials. For 
the special case where a^, a 2 and a^ are equal, equation (39) reduces 
to equation (29). 

Other approximations could be used to develop the difference 
equations at the interface. For example. 


9 

l l 9x 


- a 


9<f) 
2 9x 


8(f? 


^( a x a 2 ) 3x 


(a l V 

2 Ax 


( 42 ) 
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In this case, the coefficients in front of the <b. , . and <]>.,- . 

terms in equation (40) would be identically equal to — ^ j. 

This approximation would average the temperature gradient at the in- 
terface. It would, in effect, give the same results had the interface 
between materials 1 and 2 been positioned midway between the grid 
points. With the interface positioned midway between the grid points, 
the difference equations to the left of the interface would be 
written as 


Ay) *i,j-l + Vl,j ~ 2 


1 + r) + 1 — (2ttA x) 2 

^y J a i 




+ ^ 


i+i.J +^ 2 *l,j+l " ° 


W) 


while to the right of the interface 


Ax 

fly. 


> ±, j-l + ^i-l , j 


- 2 




AyJ 


(2-irAx)' 






A 2 


+ h + i,j + (If; *i.3+i = 0 


(44) 


If equation (43) were multiplied by — , equation (44) multiplied by 

a 2 r 1 / a i +a 2'\ 

— and the resulting equations averaged, then — 1 - J would be the 

coefficient on the $ terms, which would be consistant with equation (42), 

With sufficient number of grid points, either approximation should 

yield approximately the same temperature distribution except vefy near to 

the interface . 
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CONCLUDING REMARKS 

A finite difference solution for heat wave propagation in a two- 
dimensional medium has been presented. A special exit impedance is 
derived which allows a finite heat field to approximate a semi-infinite 
medium. The numerical theory Is shown to be in good agreement with 
the corresponding exact analytical results for a semi-infinite one- 
dimensional medium. 

The finite difference formulation is flexible and should be a 
useful tool in the solution of complex heat conduction problems . 
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Figure 4. - Analytical and numerical temperature 
variations for one-dimensional heat wave prop- 
agation into a semi-infinite medium. 
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